#pragma once

#include "vector.h"

namespace g2 {

template <typename T>
class basic_triangle;

typedef basic_triangle<double> triangle;
typedef basic_triangle<float> trianglef;
typedef basic_triangle<long double> trianglel;

template <typename T>
class basic_triangle {
public:
	basic_vector<T> v[3]; ///< vertices

public:
	basic_triangle()
		{ }
	basic_triangle(const basic_vector<T> & v0, const basic_vector<T> & v1, const basic_vector<T> & v2)
		{ v[0] = v0; v[1] = v1; v[2] = v2; }
	template <class U> 
	basic_triangle(const basic_triangle<U> & t)
		{ v[0] = T(t.v[0]); v[1] = T(t.v[1]); v[2] = T(t.v[2]); }
	basic_vector<T> & operator [](size_t n)
		{ return v[n]; }
	const basic_vector<T> & operator [](size_t n) const
		{ return v[n]; }
	bool operator == (const basic_triangle<T> & t) const
		{ return v[0] == t.v[0] && v[1] == t.v[1] && v[2] == t.v[2]; }
	bool operator != (const basic_triangle<T> & t) const
		{ return !(*this == t); }
	bool operator < (const basic_triangle<T> & t) const ///<lexicographical order
		{ return v[0] < t.v[0] || v[0] == t.v[0] && 
			(v[1] < t.v[1] || v[1] == t.v[1] && v[2] < t.v[2]); }
	basic_triangle<T> & operator = (const basic_triangle<T> & t) 
		{ v[0] = t.v[0]; v[1] = t.v[1]; v[2] = t.v[2]; return * this; }

	///signed doubled value of basic_triangle's area, negative value is returned for not ccw order of vertices
	T double_signed_area() const
		{ return outer_product(v[2]-v[1], v[0]-v[1]); }
	T area() const ///< positive areas of the triangle, slower than double_signed_area()
		{ return T(0.5) * fabs(double_signed_area()); }
	///return true for counter-clockwise following of vertices
	bool ccw() const
		{ return double_signed_area() > 0; }
	basic_triangle<T> & make_ccw()
		{ if (!ccw()) std::swap(v[0], v[1]); return * this; }
	T perimeter() const
		{ return absolute(v[0]-v[1]) + absolute(v[1]-v[2]) + absolute(v[2]-v[0]); }

public:	//complex operatrion
	///checks whether the point is inside the triangle, the triangle must be ccw
	bool inside(const basic_vector<T> & p) const;
	bool inside_or_on_bound(const basic_vector<T> & p) const;
	///return signed height (positive for ccw), dropped from node #n
	T height(int n) const;
	T squared_height(int n) const; ///< sligthly faster than height(...)
	///return projection of node #n on opposite edge
	basic_vector<T> projection(int n) const;
	///retuns height, dropped from node #n, if it falls inside opposite side of the triangle,
	///otherwise returns the smaller distance to one of other vertices;
	///optionally the functions return the nearest point from the opposite edge (either projection on it or one of its ends)
	T distance_to_other_size(int n, basic_vector<T> * nearest = 0) const; ///< slower that squared_... version of the method
	T squared_distance_to_other_size(int n, basic_vector<T> * nearest = 0) const;

	///\return on demand radius and/or center of the inscribed circle
	void incircle(T * radius, basic_vector<T> * center) const;
	///\return radius of the inscribed circle, its negative for not ccw order of vertices
	T incircle_radius() const;
	///\return center of the inscribed circle
	basic_vector<T> incircle_center() const;

	///\return on demand radius and/or center of the inscribed half-circle with center on edge #e (opposite to node #e)
	void inhalfcircle(int e, T * radius, basic_vector<T> * center) const;
	///\return radius of the inscribed half-circle, its negative for not ccw order of vertices
	T inhalfcircle_radius(int e) const;
	///\return center of the inscribed half-circle
	basic_vector<T> inhalfcircle_center(int e) const;

	///\return on demand radius and/or center of the circumscribing circle
	void circumcircle(T * radius, basic_vector<T> * center, T * squared_radius = 0) const;
	///\return radius of the circumscribing circle
	T circumcircle_radius() const;
	T circumcircle_squared_radius() const;
	///\return center of the circumscribing circle
	basic_vector<T> circumcircle_center() const;

	///\return on demand radius and/or center of the special circumscribing circle with center somewhere on edge #e (opposite to node #e),
	///this circle lies totally inside the normal circumcircle, touching it at one point
	void circumhalfcircle(int e, T * radius, basic_vector<T> * center) const;
	///\return radius of the circumscribing half-circle
	T circumhalfcircle_radius(int e) const;
	///\return center of the circumscribing half-circle
	basic_vector<T> circumhalfcircle_center(int e) const;

public: //interpolation
	///computes relative areas of the triangles formed with the given point:
	///s0 + s1 + s2 = 1
	void interpolation_coordinates(basic_vector<T> const & p, T & s0, T & s1, T & s2) const;
	///given values at the vertices of the triangle reconstructs value at the given point
	///using linear interpolation if the point is inside the triangle or linear extrapolation otherwise
	template <typename U>
	U linear(basic_vector<T> const & p, U const & u0, U const & u1, U const & u2) const;
	///given values at the vertices (u0,u1,u2) and middles of edges (u01, u12, u21) of the triangle,
	///reconstructs value at the given point using quadratic interpolation 
	///if the point is inside the triangle or quadratic extrapolation otherwise
	template <typename U>
	U quadratic(basic_vector<T> const & p, U const & u0, U const & u1, U const & u2,
		U const & u01, U const & u12, U const & u20) const;
	template <typename U>
	static U quadratic_impl(T s0, T s1, T s2, U const & u0, U const & u1, U const & u2,
		U const & u01, U const & u12, U const & u20);
	///given the same kind of values as for quadratic(...),
	///reconstructs value at the given point using piecewise linear interpolation
	///(linear in each of four subtriangles)
	template <typename U>
	U piecewise_linear(basic_vector<T> const & p, U const & u0, U const & u1, U const & u2,
		U const & u01, U const & u12, U const & u20) const;
	///quadratic interpolation with limited peaks, 
	///that does not increase total variation of reconstructed function over piecewise constant interpolation
	template <typename U>
	U tv0_quadratic(basic_vector<T> const & p, U const & u0, U const & u1, U const & u2,
		U const & u01, U const & u12, U const & u20) const;
	///combination of quadratic and piecewise linear interpolation;
	///predicate recives derivatives along direction in two points and must return true if qudratic function has extremum between that points
	template <typename U, typename P>
	U monotone_quadratic(basic_vector<T> const & p, 
		U const & u0, U const & u1, U const & u2, U const & u01, U const & u12, U const & u20,
		P pred) const;
	///this version of the function can work with scalars in nodes only
	T monotone_quadratic(basic_vector<T> const & p, T u0, T u1, T u2, T u01, T u12, T u20) const;
	///given values at the vertices (u0,u1,u2), on the sides (u001, u110, u112, u221, u220, u002)
	///and at the center of mass (u3) of the triangle, reconstructs value at the given point using cubic interpolation 
	///if the point is inside the triangle or cubic extrapolation otherwise
	template <typename U>
	U cubic(basic_vector<T> const & p, U const & u0, U const & u1, U const & u2,
		U const & u001, U const & u110, U const & u112, U const & u221, 
		U const & u220, U const & u002, U const & u3) const;
	template <typename U>
	static U cubic_impl(T s0, T s1, T s2, U const & u0, U const & u1, U const & u2,
		U const & u001, U const & u110, U const & u112, U const & u221, 
		U const & u220, U const & u002, U const & u3);
	///cubic interpolation with limited peaks, 
	///that practically does not increase total variation of reconstructed function over piecewise constant interpolation
	template <typename U>
	U tv0_cubic(basic_vector<T> const & p, U const & u0, U const & u1, U const & u2,
		U const & u001, U const & u110, U const & u112, U const & u221, 
		U const & u220, U const & u002, U const & u3) const;
	///given the same kind of values as for cubic(...),
	///reconstructs value at the given point using piecewise linear interpolation
	///(linear in each of nine subtriangles)
	template <typename U>
	U piecewise_linear(basic_vector<T> const & p, U const & u0, U const & u1, U const & u2,
		U const & u001, U const & u110, U const & u112, U const & u221, 
		U const & u220, U const & u002, U const & u3) const;
};

//implementation

template <typename T>
inline bool basic_triangle<T>::inside(const basic_vector<T> & p) const
{
	assert (ccw());
	return 
		outer_product(v[1] - v[0], p - v[0]) > 0 &&
		outer_product(v[2] - v[1], p - v[1]) > 0 &&
		outer_product(v[0] - v[2], p - v[2]) > 0;
}

template <typename T>
inline bool basic_triangle<T>::inside_or_on_bound(const basic_vector<T> & p) const
{
	assert (ccw());
	return 
		outer_product(v[1] - v[0], p - v[0]) >= 0 &&
		outer_product(v[2] - v[1], p - v[1]) >= 0 &&
		outer_product(v[0] - v[2], p - v[2]) >= 0;
}

template <typename T>
inline T basic_triangle<T>::height(int n) const
{ 
	int n1 = (n+1)%3;
	int n2 = (n+2)%3;
	basic_vector<T> AB = v[n2] - v[n1];
	return outer_product(AB, v[n] - v[n1]) / absolute(AB);
}

template <typename T>
inline T basic_triangle<T>::squared_height(int n) const
{ 
	int n1 = (n+1)%3;
	int n2 = (n+2)%3;
	basic_vector<T> AB = v[n2] - v[n1];
	return sqr(outer_product(AB, v[n] - v[n1])) / sqr(AB);
}

template <typename T>
inline basic_vector<T> basic_triangle<T>::projection(int n) const
{
	int n1 = (n+1)%3;
	int n2 = (n+2)%3;
	basic_vector<T> AB = v[n2] - v[n1];
	return v[n1] + AB * (scalar_product(AB, v[n] - v[n1]) / sqr(AB));
}

template <typename T>
inline T basic_triangle<T>::distance_to_other_size(int n, basic_vector<T> * nearest) const
	{ return sqrt(squared_distance_to_other_size(n, nearest)); }

template <typename T>
inline T basic_triangle<T>::squared_distance_to_other_size(int n, basic_vector<T> * nearest) const
{
	int n1 = (n+1)%3;
	int n2 = (n+2)%3;
	//old tested solution:
	//if (scalar_product(v[n1] - v[n], v[n2] - v[n1]) * scalar_product(v[n2] - v[n], v[n2] - v[n1]) < 0)
	//	return squared_height(n);
	//else
	//	return std::min(sqr(v[n1] - v[n]), sqr(v[n2] - v[n]));
	//new one:
	if (scalar_product(v[n] - v[n1], v[n2] - v[n1]) < 0) {
		if (nearest) *nearest = v[n1];
		return sqr(v[n] - v[n1]);
	}
	if (scalar_product(v[n] - v[n2], v[n1] - v[n2]) < 0) {
		if (nearest) *nearest = v[n2];
		return sqr(v[n] - v[n2]);
	}
	if (nearest) {
		*nearest = projection(0);
		return sqr(v[n] - *nearest);
	}
	return squared_height(n);
}

template <typename T>
inline void basic_triangle<T>::incircle(T * radius, basic_vector<T> * center) const
{
	if (!radius && !center)
		return;
	T d01 = absolute(v[0] - v[1]);
	T d02 = absolute(v[0] - v[2]);
	T d12 = absolute(v[1] - v[2]);
	T h = double_signed_area() / (d01 + d02 + d12);
	if (radius)
		*radius = h;
	if (!center)
		return;
	basic_vector<T> t = (v[1] - v[0]) / d01;
	basic_vector<T> n(-t.y, t.x);
	*center = v[0] + T(0.5)*(d01 + d02 - d12)*t + h*n;
}

template <typename T>
inline T basic_triangle<T>::incircle_radius() const
{
	T radius;
	incircle(&radius, 0);
	return radius;
}

template <typename T>
inline basic_vector<T> basic_triangle<T>::incircle_center() const
{
	basic_vector<T> center;
	incircle(0, &center);
	return center;
}

template <typename T>
inline void basic_triangle<T>::inhalfcircle(int e, T * radius, basic_vector<T> * center) const
{
	if (!radius && !center)
		return;
	int e1 = (e+1)%3;
	int e2 = (e+2)%3;
	basic_vector<T> AB = v[e2] - v[e1];
	basic_vector<T> AC = v[e] - v[e1];
	basic_vector<T> BC = v[e] - v[e2];
	T dAC = absolute(AC);
	T dBC = absolute(BC);
	T r = 1.0 / (dAC + dBC);
	if (radius)
		*radius = outer_product(AB, AC) * r;
	if (center)
		*center = v[e1] + (dAC * r) * AB;
}

template <typename T>
inline T basic_triangle<T>::inhalfcircle_radius(int e) const
{
	T radius;
	inhalfcircle(e, &radius, 0);
	return radius;
}

template <typename T>
inline basic_vector<T> basic_triangle<T>::inhalfcircle_center(int e) const
{
	basic_vector<T> center;
	inhalfcircle(e, 0, &center);
	return center;
}


template <typename T>
inline void basic_triangle<T>::circumcircle(T * radius, basic_vector<T> * center, T * squared_radius) const
{
	T b1 = sqr(v[0]) - sqr(v[1]);
	T b2 = sqr(v[0]) - sqr(v[2]);
	basic_vector<T> d1 = T(2) * (v[0] - v[1]);
	basic_vector<T> d2 = T(2) * (v[0] - v[2]);
	basic_vector<T> ort_d1(d1.y, -d1.x);
	basic_vector<T> ort_d2(d2.y, -d2.x);
	T u = scalar_product(d1, ort_d2);
	T w = scalar_product(ort_d1, d2);
	basic_vector<T> c;
	if (u != 0 && w != 0) {
		basic_vector<T> e1 = 1 / u * ort_d2;
		basic_vector<T> e2 = 1 / w * ort_d1;
		c = b1 * e1 + b2 * e2;
	} else 
		c = v[0]; ///< triangle is singular
	if (center)
		*center = c;
	T RR = sqr(v[0] - c);
	if (radius)
		*radius = sqrt(RR);
	if (squared_radius)
		*squared_radius = RR;
}

template <typename T>
inline T basic_triangle<T>::circumcircle_radius() const
{
	T radius;
	circumcircle(&radius, 0, 0);
	return radius;
}

template <typename T>
inline T basic_triangle<T>::circumcircle_squared_radius() const
{
	T squared_radius;
	circumcircle(0, 0, &squared_radius);
	return squared_radius;
}

template <typename T>
inline basic_vector<T> basic_triangle<T>::circumcircle_center() const
{
	basic_vector<T> center;
	circumcircle(0, &center, 0);
	return center;
}

template <typename T>
inline void basic_triangle<T>::circumhalfcircle(int e, T * radius, basic_vector<T> * center) const
{
	T R;
	basic_vector<T> C;
	circumcircle(&R, &C);
	int e1 = (e+1)%3;
	int e2 = (e+2)%3;
	basic_vector<T> d = v[e2] - v[e1];
	C -= v[e1];
	T sc = scalar_product(C,d);
	T rec = T(1) / sqr(d);
	if (center)
		*center = v[e1] + sc * rec * d;
	if (radius) {
		if (outer_product(d, C) < 0)
			//center of standard circumcircle is outside the body
			*radius = R - sqrt(sqr(C) - sqr(sc) * rec);
		else
			*radius = sqrt(sqr(R) - sqr(C) + sqr(sc) * rec);
	}
}

template <typename T>
inline T basic_triangle<T>::circumhalfcircle_radius(int e) const
{
	T radius;
	circumhalfcircle(e, &radius, 0);
	return radius;
}

template <typename T>
inline basic_vector<T> basic_triangle<T>::circumhalfcircle_center(int e) const
{
	basic_vector<T> center;
	circumhalfcircle(e, 0, &center);
	return center;
}

// interpolation

template <typename T>
inline void basic_triangle<T>::interpolation_coordinates(basic_vector<T> const & p, T & s0, T & s1, T & s2) const
{
	s2 = outer_product(v[1] - v[0], p - v[0]);
	s0 = outer_product(v[2] - v[1], p - v[1]);
	s1 = outer_product(v[0] - v[2], p - v[2]);
	T rS = 1 / (s0 + s1 + s2);
	s0 *= rS;
	s1 *= rS;
	s2 *= rS;
}

template <typename T>
template <typename U>
inline U basic_triangle<T>::linear(basic_vector<T> const & p, U const & u0, U const & u1, U const & u2) const
{
	T s0, s1, s2;
	interpolation_coordinates(p, s0, s1, s2);
	return u0 * s0 + u1 * s1 + u2 * s2;
}

template <typename T>
template <typename U>
inline U basic_triangle<T>::quadratic(basic_vector<T> const & p, U const & u0, U const & u1, U const & u2,
	U const & u01, U const & u12, U const & u20) const
{
	T s0, s1, s2;
	interpolation_coordinates(p, s0, s1, s2);
	return quadratic_impl(s0, s1, s2, u0, u1, u2, u01, u12, u20);
}

template <typename T>
template <typename U>
inline U basic_triangle<T>::quadratic_impl(T s0, T s1, T s2, U const & u0, U const & u1, U const & u2,
	U const & u01, U const & u12, U const & u20)
{
	return 
		s0 * (s0 - s1 - s2) * u0 + 
		s1 * (s1 - s0 - s2) * u1 + 
		s2 * (s2 - s0 - s1) * u2 +
		4 * s0 * s1 * u01 +
		4 * s1 * s2 * u12 +
		4 * s2 * s0 * u20;
}

template <typename T>
template <typename U>
inline U basic_triangle<T>::piecewise_linear(basic_vector<T> const & p, U const & u0, U const & u1, U const & u2,
	U const & u01, U const & u12, U const & u20) const
{
	T s0, s1, s2;
	interpolation_coordinates(p, s0, s1, s2);
	s0 *= 2;
	s1 *= 2;
	s2 *= 2;
	if (s0 >= 1) 
		return (s0 - 1) * u0 + s1 * u01 + s2 * u20;
	if (s1 >= 1) 
		return (s1 - 1) * u1 + s0 * u01 + s2 * u12;
	if (s2 >= 1) 
		return (s2 - 1) * u2 + s0 * u20 + s1 * u12;
	return (1 - s2) * u01 + (1 - s0) * u12 + (1 - s1) * u20;
}

template <typename T>
template <typename U>
U basic_triangle<T>::tv0_quadratic(basic_vector<T> const & p, U const & u0, U const & u1, U const & u2,
	U const & u01, U const & u12, U const & u20) const
{
	T s0, s1, s2;
	interpolation_coordinates(p, s0, s1, s2);
	U q = quadratic_impl(s0, s1, s2, u0, u1, u2, u01, u12, u20);
	U mn, mx;
	if (2*s0 >= 1) {
		mn = std::min(u0, std::min(u01, u20));
		mx = std::max(u0, std::max(u01, u20));
	}
	else if (2*s1 >= 1) {
		mn = std::min(u1, std::min(u01, u12));
		mx = std::max(u1, std::max(u01, u12));
	}
	else if (2*s2 >= 1) {
		mn = std::min(u2, std::min(u20, u12));
		mx = std::max(u2, std::max(u20, u12));
	}
	else {
		mn = std::min(u01, std::min(u12, u20));
		mx = std::max(u01, std::max(u12, u20));
	}
	return std::max(std::min(q, mx), mn);
}

template <typename T>
template <typename U, typename P>
U basic_triangle<T>::monotone_quadratic(basic_vector<T> const & p, 
	U const & u0, U const & u1, U const & u2, U const & u01, U const & u12, U const & u20,
	P pred) const
{
	//basic_vector<T> k0 = (v[2] - v[1]).orthogonal();
	//basic_vector<T> k1 = (v[0] - v[2]).orthogonal();
	//basic_vector<T> k2 = (v[1] - v[0]).orthogonal();

	T s0, s1, s2;
	interpolation_coordinates(p, s0, s1, s2);
	s0 *= 2;
	s1 *= 2;
	s2 *= 2;

	//gradients of the interpolated function at the vertices of the triangle
	//basic_vector<T> U0 = 3 * u0 * k0 + (4 * u01 - u1) * k1 + (4 * u20 - u2) * k2;
	//basic_vector<T> U1 = 3 * u1 * k1 + (4 * u01 - u0) * k0 + (4 * u12 - u2) * k2;
	//basic_vector<T> U2 = 3 * u2 * k2 + (4 * u20 - u0) * k0 + (4 * u12 - u1) * k1;

	//scalar products of the gradient in nodes on the edge directions
	U sc_U0_l0 = 4 * (u20 - u01) - (u2 - u1);
	U sc_U1_l1 = 4 * (u01 - u12) - (u0 - u2);
	U sc_U2_l2 = 4 * (u12 - u20) - (u1 - u0);
	U sc_U0_l1 =  3 * u0 - 4 * u20 + u2;
	U sc_U0_l2 = -3 * u0 + 4 * u01 - u1;
	U sc_U1_l2 =  3 * u1 - 4 * u01 + u0;
	U sc_U1_l0 = -3 * u1 + 4 * u12 - u2;
	U sc_U2_l0 =  3 * u2 - 4 * u12 + u1;
	U sc_U2_l1 = -3 * u2 + 4 * u20 - u0;

	if (s0 >= 1) {
		// 0 - 01 - 20 subtriangle
		U u_0_01 = //(outer_product(k2, U0) * outer_product(k2, U0 + U1) < 0) ?
			pred(sc_U0_l2, sc_U0_l2 + sc_U1_l2) ?
			T(0.5) * (u0 + u01) : //linear
			T(0.375) * u0 + T(0.75) * u01 - T(0.125) * u1; //quadratic
		U u_01_20 = //(outer_product(k0, U0 + U1) * outer_product(k0, U0 + U2) < 0) ?
			pred(sc_U0_l0 + sc_U1_l0, sc_U0_l0 + sc_U2_l0) ?
			T(0.5) * (u01 + u20) : //linear
			T(0.5) * (u01 + u20) + T(0.25) * u12 - T(0.125) * (u1 + u2); //quadratic
		U u_20_0 = //(outer_product(k1, U0) * outer_product(k1, U0 + U2) < 0) ?
			pred(sc_U0_l1, sc_U0_l1 + sc_U2_l1) ?
			T(0.5) * (u0 + u20) : //linear
			T(0.375) * u0 + T(0.75) * u20 - T(0.125) * u2; //quadratic
		return quadratic_impl(s0 - 1, s1, s2, u0, u01, u20, u_0_01, u_01_20, u_20_0);
	}
	else if (s1 >= 1) {
		// 1 - 12 - 01 subtriangle
		U u_1_12 = //(outer_product(k0, U1) * outer_product(k0, U1 + U2) < 0) ?
			pred(sc_U1_l0, sc_U1_l0 + sc_U2_l0) ?
			T(0.5) * (u1 + u12) : //linear
			T(0.375) * u1 + T(0.75) * u12 - T(0.125) * u2; //quadratic
		U u_12_01 = //(outer_product(k1, U1 + U2) * outer_product(k1, U1 + U0) < 0) ?
			pred(sc_U1_l1 + sc_U2_l1, sc_U1_l1 + sc_U0_l1) ?
			T(0.5) * (u12 + u01) : //linear
			T(0.5) * (u12 + u01) + T(0.25) * u20 - T(0.125) * (u0 + u2); //quadratic
		U u_01_1 = //outer_product(k2, U1) * (outer_product(k2, U1 + U0) < 0) ?
			pred(sc_U1_l2, sc_U1_l2 + sc_U0_l2) ?
			T(0.5) * (u1 + u01) : //linear
			T(0.375) * u1 + T(0.75) * u01 - T(0.125) * u0; //quadratic
		return quadratic_impl(s1 - 1, s2, s0, u1, u12, u01, u_1_12, u_12_01, u_01_1);
	}
	else if (s2 >= 1) {
		// 2 - 20 - 12 subtriangle
		U u_2_20 = //(outer_product(k1, U2) * outer_product(k1, U2 + U0) < 0) ?
			pred(sc_U2_l1, sc_U2_l1 + sc_U0_l1) ?
			T(0.5) * (u2 + u20) : //linear
			T(0.375) * u2 + T(0.75) * u20 - T(0.125) * u0; //quadratic
		U u_20_12 = //(outer_product(k2, U1 + U2) * outer_product(k2, U2 + U0) < 0) ?
			pred(sc_U1_l2 + sc_U2_l2, sc_U2_l2 + sc_U0_l2) ?
			T(0.5) * (u20 + u12) : //linear
			T(0.5) * (u20 + u12) + T(0.25) * u01 - T(0.125) * (u0 + u1); //quadratic
		U u_12_2 = //(outer_product(k0, U2) * outer_product(k0, U2 + U1) < 0) ?
			pred(sc_U2_l0, sc_U2_l0 + sc_U1_l0) ?
			T(0.5) * (u12 + u2) : //linear
			T(0.375) * u2 + T(0.75) * u12 - T(0.125) * u1; //quadratic
		return quadratic_impl(s2 - 1, s0, s1, u2, u20, u12, u_2_20, u_20_12, u_12_2);
	}
	else {
		// 01 - 12 - 20 subtriangle
		U u_01_20 = //(outer_product(k0, U0 + U1) * outer_product(k0, U0 + U2) < 0) ?
			pred(sc_U0_l0 + sc_U1_l0, sc_U0_l0 + sc_U2_l0) ?
			T(0.5) * (u01 + u20) : //linear
			T(0.5) * (u01 + u20) + T(0.25) * u12 - T(0.125) * (u1 + u2); //quadratic
		U u_12_01 = //(outer_product(k1, U1 + U2) * outer_product(k1, U1 + U0) < 0) ?
			pred(sc_U1_l1 + sc_U2_l1, sc_U1_l1 + sc_U0_l1) ?
			T(0.5) * (u12 + u01) : //linear
			T(0.5) * (u12 + u01) + T(0.25) * u20 - T(0.125) * (u0 + u2); //quadratic
		U u_20_12 = //(outer_product(k2, U1 + U2) * outer_product(k2, U2 + U0) < 0) ?
			pred(sc_U1_l2 + sc_U2_l2, sc_U2_l2 + sc_U0_l2) ?
			T(0.5) * (u20 + u12) : //linear
			T(0.5) * (u20 + u12) + T(0.25) * u01 - T(0.125) * (u0 + u1); //quadratic
		return quadratic_impl(1 - s0, 1 - s1, 1 - s2, u12, u20, u01, u_20_12, u_01_20, u_12_01);
	}
}

template <typename T>
struct scalar_derivative_zero_test {
	bool operator() (T d0, T d1) const
		{ return d0 * d1 < 0; }
};

template <typename T>
inline T basic_triangle<T>::monotone_quadratic(basic_vector<T> const & p, T u0, T u1, T u2,
	T u01, T u12, T u20) const
	{ return monotone_quadratic(p, u0, u1, u2, u01, u12, u20, scalar_derivative_zero_test<T>()); }

template <typename T>
template <typename U>
inline U basic_triangle<T>::cubic_impl(T s0, T s1, T s2, U const & u0, U const & u1, U const & u2,
	U const & u001, U const & u110, U const & u112, U const & u221, U const & u220, U const & u002, 
	U const & u3)
{
	T w0, w1, w2, w001, w110, w112, w221, w220, w002, w3;
	w0 = T(0.5) * s0 * (3 * s0 - 1) * (3 * s0 - 2);
	w1 = T(0.5) * s1 * (3 * s1 - 1) * (3 * s1 - 2);
	w2 = T(0.5) * s2 * (3 * s2 - 1) * (3 * s2 - 2);
	w001 = T(4.5) * s0 * s1 * (3*s0 - 1);
	w110 = T(4.5) * s0 * s1 * (3*s1 - 1);
	w112 = T(4.5) * s1 * s2 * (3*s1 - 1);
	w221 = T(4.5) * s1 * s2 * (3*s2 - 1);
	w220 = T(4.5) * s2 * s0 * (3*s2 - 1);
	w002 = T(4.5) * s2 * s0 * (3*s0 - 1);
	w3 = 27 * s1 * s2 * s0;
	return w0*u0 + w1*u1 + w2*u2 + w001*u001 + w110*u110 + w112*u112 +
		w221*u221 + w220*u220 + w002*u002 + w3*u3;
}

template <typename T>
template <typename U>
inline U basic_triangle<T>::cubic(basic_vector<T> const & p, U const & u0, U const & u1, U const & u2,
	U const & u001, U const & u110, U const & u112, U const & u221, U const & u220, U const & u002, 
	U const & u3) const
{
	T s0, s1, s2;
	interpolation_coordinates(p, s0, s1, s2);
	return cubic_impl(s0, s1, s2, u0, u1, u2,
		u001, u110, u112, u221, u220, u002, u3);
}

template <typename T>
template <typename U>
inline U basic_triangle<T>::piecewise_linear(basic_vector<T> const & p, U const & u0, U const & u1, U const & u2,
	U const & u001, U const & u110, U const & u112, U const & u221, 
	U const & u220, U const & u002, U const & u3) const
{
	T s0, s1, s2;
	interpolation_coordinates(p, s0, s1, s2);
	s0 *= 3;
	s1 *= 3;
	s2 *= 3;
	//0 - 001 - 002
	if (s0 >= 2)
		return (s0 - 2) * u0 + s1 * u001 + s2 * u002;
	//1 - 112 - 110
	if (s1 >= 2) 
		return (s1 - 2) * u1 + s0 * u110 + s2 * u112;
	//2 - 220 - 221
	if (s2 >= 2) 
		return (s2 - 2) * u2 + s0 * u220 + s1 * u221;
	//112 - 221 - 3
	if (s1 >= 1 && s2 >= 1) 
		return s0 * u3 + (s2 - 1) * u221 + (s1 - 1) * u112;
	//002 - 220 - 3
	if (s2 >= 1 && s0 >= 1) 
		return (s0 - 1) * u002 + s1 * u3 + (s2 - 1) * u220;
	//110 - 001 - 3
	if (s1 >= 1 && s0 >= 1)
		return (s1 - 1) * u001 + (s0 - 1) * u110 + s2 * u3;
	//221 - 220 - 3
	if (s0 <= 1 && s1 <= 1)
		return (1 - s1) * u220 + (1 - s0) * u221 + (2 - s2) * u3;
	//110 - 112 - 3
	if (s0 <= 1 && s2 <= 1)
		return (1 - s2) * u110 + (2 - s1) * u3 + (1 - s0) * u112;
	//002 - 001 - 3
	return (2 - s0) * u3  + (1 - s2) * u001 + (1 - s1) * u002;
}

template <typename T>
template <typename U>
U basic_triangle<T>::tv0_cubic(basic_vector<T> const & p, U const & u0, U const & u1, U const & u2,
	U const & u001, U const & u110, U const & u112, U const & u221, U const & u220, U const & u002, 
	U const & u3) const
{
	T s0, s1, s2;
	interpolation_coordinates(p, s0, s1, s2);
	U q = cubic_impl(s0, s1, s2, u0, u1, u2,
		u001, u110, u112, u221, u220, u002, u3);
	U mn, mx;
	if (3*s0 >= 2) {
		mn = std::min(u0, std::min(u001, u002));
		mx = std::max(u0, std::max(u001, u002));
	}
	else if (3*s1 >= 2) {
		mn = std::min(u1, std::min(u110, u112));
		mx = std::max(u1, std::max(u110, u112));
	}
	else if (3*s2 >= 2) {
		mn = std::min(u2, std::min(u220, u221));
		mx = std::max(u2, std::max(u220, u221));
	}
	else if (3*s1 >= 1 && 3*s2 >= 1) {
		mn = std::min(u3, std::min(u112, u221));
		mx = std::max(u3, std::max(u112, u221));
	}
	else if (3*s2 >= 1 && 3*s0 >= 1) {
		mn = std::min(u3, std::min(u220, u002));
		mx = std::max(u3, std::max(u220, u002));
	}
	else if (3*s1 >= 1 && 3*s0 >= 1) {
		mn = std::min(u3, std::min(u001, u110));
		mx = std::max(u3, std::max(u001, u110));
	}
	else if (3*s0 <= 1 && 3*s1 <= 1) {
		mn = std::min(u3, std::min(u220, u221));
		mx = std::max(u3, std::max(u220, u221));
	}
	else if (3*s0 <= 1 && 3*s2 <= 1) {
		mn = std::min(u3, std::min(u110, u112));
		mx = std::max(u3, std::max(u110, u112));
	}
	else {
		mn = std::min(u3, std::min(u002, u001));
		mx = std::max(u3, std::max(u002, u001));
	}
	return std::max(std::min(q, mx), mn);
}

} // g2
